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Abstract 


An efficient direct and numerical method has been proposed to approx- 
imate a solution of time-delay fractional optimal control problems. First, a 
class of discrete orthogonal polynomials, called Hahn polynomials, has been 
introduced and their properties are investigated. These properties are em- 
ployed to derive a general formulation of their operational matrix of fractional 
integration, in the Riemann—Liouville sense. Then, the fractional derivative 
of the state function in the dynamic constraint of time-delay fractional op- 
timal control problems is approximated by the Hahn polynomials with un- 
known coefficients. The operational matrix of fractional integration together 
with the dynamical constraints is used to approximate the control function 
directly as a function of the state function. Finally, these approximations 
were put in the performance index and necessary conditions for optimality 
transform the under consideration time-delay fractional optimal control prob- 
lems into an algebraic system. Some illustrative examples are given and the 
obtained numerical results are compared with those previously published in 
the literature. 
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1 Introduction 


It has been found that integer-order calculus is not an appropriate tool for 
modeling complex systems in science and engineering. Fractional calculus, 
as an extension of derivatives and integrals to noninteger orders, has been 
recently used to model many fundamental problems; see [32,38]. Recently, 
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fractional differential equations have found many applications in various fields 
of engineering and physics such as colored noise, signal processing, electro- 
magnetism, electrochemistry, dynamic of viscoelastic materials, continuum 
and statistical mechanics, solid mechanics, and fluid-dynamic traffic model; 
see [3,4,9, 12, 40]. 


Contrary to ordinary differential equations (ODEs), delay differential 
equations (DDEs) contain derivatives that depend on the solution at previ- 
ous times, thus making the mathematical model closer to the real-world phe- 
nomenon; see [26,27]. Recently, many researches have been focused on DDEs 
and their numerical solution. For example, the finite difference method [22], 
Adams methods [41], Adomian decomposition method [13], variational itera- 
tion method [45], hybrid functions [29], spectral methods [19, 39, 44], and 
wavelet methods [35,36] have been employed to approximate solution of 
DDEs. 


Over the last decades, numerical methods based on orthogonal functions 
have been frequently used to approximate solution of differential and integral 
equations [7,28,31,34,39]. The main characteristic of these methods is that 
they reduce under consideration problems to those of solving a system of al- 
gebraic equations, thus greatly simplifying the problems. Depending on their 
structure, the orthogonal functions are classified into three families: the first 
class consists of sets of piecewise constant orthogonal functions such as the 
Walsh functions and block pulse functions. The second class consists of sets 
of sine-cosine functions, and the last class is the most widely used orthog- 
onal polynomials, such as Laguerre, Legendre, and Chebyshev polynomials; 
see [28,31]. It is well known that orthogonal polynomials that are solutions of 
singular Sturm—Liouville problems allow approximation of smooth functions, 
where truncation error approaches zero faster than any negative power of the 
number of basic functions used in the approximation. This phenomenon is 
usually referred to as spectral accuracy; see [8,30]. Moreover, according to 
the defined inner product in the solution space, orthogonal polynomials are 
classified into two main classes: continuous and discrete. For continuous or- 
thogonal polynomials such as Legendre, Chebyshev, Hermite, and Laguerre, 
one has to evaluate an integral in the inner product, whereas discrete or- 
thogonal polynomials come with a discrete scalar product and hence the 
integral becomes a sum. Although continuous orthogonal polynomials have 
been more frequently used to approximate solution of functional equations, 
there are some advantages of using discrete orthogonal polynomials. For ex- 
ample, by using discrete orthogonal polynomials, the Fourier coefficients can 
be calculated with the aid of a summation and the obtained coefficients are 
exact. Consequently, comparison to continuous cases, implementation of dis- 
crete orthogonal polynomials is more efficient and less complex; see [16, 43]. 
Moreover, there is a close connection between stochastic processes and dis- 
crete orthogonal polynomials that motivated many researchers to consider 
them to solve stochastic differential equations; see [2,43]. 
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The optimal control of a system is one of the most practical subjects in 
science and engineering [5, 6, 17, 21, 25,33]. The optimal control problems 
involve minimization of a performance index subject to a dynamical system. 
As generalizations of the classical optimal control, FOCPs are problems in 
which fractional derivatives or integrals are used in the performance index 
or constraints. Similar to the other types of fractional functional equations, 
most of FOCPs do not have exact and analytic solutions [7,34,35,37]. Fur- 
thermore, a time-delay fractional optimal control problem (TDFOCP) is a 
fractional optimal control problem in which the dynamic of system contains 
some time-delay equations. The control of time-delay systems frequently in 
electronic, age structure, biological, chemical, electronic, and transportation 
systems [14, 20, 23,41]. Due to their variety of applications in the realistic 
models of phenomena, the control of time-delay systems has been investi- 
gated by many engineers and scientists. Also, considerable attention has 
been focused on the approximate and numerical solution of them [7,34,35,37]. 
Generally, there are two main approaches to the approximate solution of TD- 
FOCPs. The first one is based on constructing a Hamiltonian system and 
then solve the arising two-point boundary value problem, which is the indirect 
approach, and the other involves facing directly the problems by discretizing 
or approximating the functions without constructing the Hamiltonian equa- 
tions. Recently, more attention has been done to direct approaches [7, 22,35]. 


The main purpose of this work is to introduce a new kind of discrete 
orthonormal polynomials. Some properties of these discrete polynomials are 
studied and a general formulation of their operational matrix of fractional 
integration, in the Riemann—Liouville sense, is derived. Make use of these 
orthogonal polynomials and their operational matrix, an efficient direct nu- 
merical method is proposed to approximate the solution of the following 
TDFOCP: 


Min J = al [r(t)a? (t) + s(t)u? (t)] dt, (1) 


subjected to 


D’ x(t) =a(t)x(t)+b(t)u(t)+c(t)xz(t—7) + e(t)hu(t—n)+ f(t), (2) 
t € [0,1], 0<v <1, (3) 


a(t) = a : me 0), (4) 


in which x(t) and u(t) are the state and control functions, respectively. Fur- 
thermore, r(t) and s(t) are positive functions and a(t), b(t), c(t), e(t), and 
f (€) are continuous functions, gi(t) and go(t) are arbitrary known functions 
defined on the intervals [—r,0] and [—7, 0], respectively. Also, r and 7 > 0 
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are given constant delays and D’a(t) is the fractional derivative of state 
function x(t) in the Caputo sense. 


This paper is organized as follows: In Section 2, some basic definitions 
and preliminary remarks on the fractional calculus are presented. The Hahn 
polynomials and their properties are introduced in Section 3. The operational 
matrix of fractional integration and product operational matrix for the Hahn 
polynomials have been derived in Section 4. In Section 5, an efficient direct 
method is proposed to solve TDFOCPs. To confirm the efficiency and accu- 
racy of the presented method, some illustrative examples are given in Section 
6. Finally, a conclusion is given in Section 7. 


2 Definitions and preliminaries 


Although, the fractional order operators enable to model a wider class of 
problems, there does not exist a unique definition of fractional derivative. 
Among the several formulations of the generalized derivative, the Riemann— 
Liouville and Caputo definitions are the most commonly used, which can be 
described as follows [32,38]. 


Definition 1. A real function f(t),t > 0, is said to belong the space 
C., & € R if there exist a real number p > yz and a function f\(t) € C[0,0o) 
such that f(t) = ¢t? f(t). Moreover, if f(™ € C,,, then f(t) is said to belong 
the space C7;, nEN. 


Definition 2. The Riemann—Liouville fractional integration of order v > 0 
of a function f € C,,, 4 > —1, is defined as 


: v-letir)dr, ov 
ee ray f oar, > 0, 
I, y=—0. 


The Riemann-Liouville fractional operator I” has the following properties: 


m (J f(t)) = 2? I" ft), ,v2 20, 


m (I fi))=I™ ft), 1,72 20 


T(A+1 4 
p= tate x pgs. 


Definition 3. The Caputo fractional derivative of order v > 0 is defined as 


A discrete orthogonal polynomials approach for fractional ... 53 


d” f(t) 


an v=ne€EN, 


D” f(t) = 


1 f © 
dr, t>0,0<n-1 
oe T, t>0,0<n—1<v<n, 


where n is integer, t > 0, and f € CY. 


For No = {0,1,2,...}, f € Cu, w,A > —-1, andn—-1<v <n, some useful 
and practical properties of the Caputo fractional operators D” are given by 
the following expressions: 


n-1 k 
Vv V 4 t 
PD ft) =f)- Dd FOO) g, #>0, 
k=0 
D°l’ f(t) = fF), 
0 for XENo and A < »y, 
Dp _ 
TA+)) py» otherwise 
T(A—v +1) : 


For more details on fractional calculus and their applications we refer the 
reader to [32,38]. 


3 Hahn polynomials and their properties 


The Hahn polynomials are a class of orthogonal polynomials, which are in- 
troduced and investigated by Wolfgang Hahn [18]. These polynomials can 
be considered as a discrete analog of the classical Jacobi polynomials [16,42]. 
This section is devoted to a brief definition of Hahn polynomials and their 
properties. For more details and some applications of the Hahn polynomials, 
one can refer the reader to [16, 18, 24, 42, 43]. 


Definition 4. For arbitrary complex numbers a;, i = 1,2,3 and b; #0, j = 
1,2, the generalized hypergeometric series 3 F > (a1, 2, 43; 01, b2; x) is defined 
as 


— (a1), (a2), (a3), 2" 
F: >b1, bo; 2) = ) 
3 2 (41, G2, G3; 1> 23 £) om (b1);,(b2), k! 


in which a can be derived as follows: 
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(a)o = 1, 
(a), =a(at+1)---(a+k-1), k>1, 


and the series terminates if one of the numbers a, is zero or a negative integer. 


Definition 5. Let a, > —1 be real numbers and let N be a positive integer 
number. The Hahn polynomial H,,(t;a,{,N) may be defined in terms of a 
hypergeometric series 


Qn (x; a, 8, N) = 3Fo( nntat+B+1, z;a+1,—N;1) 


= les Dee , n=0,1,...,N. 


Remark 1. Hereafter, the simple notation Q,,(x) is used to denote the nth 
Hahn polynomial Q,, (x; a, 8, N). 


Orthogonality: 
The set of Hahn polynomials {Q,(x) : n = 0,1,...,N} are orthogonal on 
the interval [0, N] with respect to the following discrete norm: 


N 


(f.9)w = df (r)gryu(r), 


r=0 
where w(t) is the weight function and is defined by 
_ fate B+N-<2 
wie) = (579) nee") 
_ T(a+a¢+ 1) (6+N-2+1) 
— T(@ +I (N-24+)l(attl(B+l) 


Moreover, the orthogonality condition for these polynomials reads as follows: 


N 
(On) On B= > OnMOmr yo) =a (i) oma 
r=0 


where dmn is the Kronecker delta and a(n) can be defined as 


(-1)"(n+a4+68+1)y4, (841), 0! 


™(")= "GatatBt(atl,(-N), Nl 


Recurrence relation: 
The set of Hahn polynomials {Q,(a) : n = 0,1,...,N} can be determined 
with the aid of the following recurrence formulas: 


HnQn41(2) = (Ln + On — x) Qn(x) = OnQn-1(2), n= 0, 1, See N- ale 


A discrete orthogonal polynomials approach for fractional ... 55 


where 


(ntat+8+1)(ntat1)(N-n) 
(Qnt+a+6+4+1)(2Qn+a+8+2) 


Ln = 


? 


and 


n(n+8)(n+at+B+N) 
Qn+a+B)Qn+a+B+1)’ 


n 


Explicit formula: 
By using the generalized hypergeometric series 3/2, the explicit analytical 
form of the Hahn polynomials can be derived as follows: 


k=0 
in which 
S-Di (=n), (nt+atB+0,5 (b, 4) 
hin = 2 Nees CN: (5) 
and 


Shifted Hahn polynomial: 

In order to use the Hahn polynomials on the interval [0,1], we define the 
shifted Hahn polynomials (SHPs) by introducing the change of variable x = 
Nt. Let the nth shifted Hahn polynomial, Q,(Nt), be denoted by H,,(t). 
Then, the set of shifted Hahn polynomials {H,,(t) : n = 0,1,...,N} are 
orthogonal on the interval [0,1] with respect to the shifted weight function 
w(t) = w (Nt) and the following discrete norm: 


tbe S1(2)o(Z) (8) o 


Also, the orthogonality condition for SHPs can be defined as follows: 


(H,(t), H(t) = 3 Hy, (=) Hin (<) w (=) =n(n)6mn, (7) 
r=0 


where 7 (n) is defined by (10). 


Function approximation: 
A function f(t), defined over the interval [0, 1], may be expanded by the SHPs 
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as 
N 


f= > aki(t) = CTS), (8) 


i=0 


where C and ®(x) are (N +1) vectors given by 
CT = [co,c1,...,en], ®(t) = [Ho(t), Hi(t),..., Hn (t)]" , (9) 


and the coefficients c; can be derived as follows: 


i=0,1,...,N. (10) 


Convergence analysis: 


Convergence analysis and spectral accuracy of the Hahn polynomials expan- 
sion are investigated thoroughly in [15,16]. The following Theorem provides 
conditions that ensure that the series expansion of a function by the discrete 
Hahn polynomials converges. 


(Qk Qk) 
the discrete Hahn polynomials converges pointwise, if the series expansion 
of the function f by the Jacobi polynomials converges pointwise and if 


3+a+maz{1,c} 
x 7 0 as n, N > 00. 


Theorem 1. The series expansion >> (fn) Q(x) of a function f by 
k=0 


Proof. The proof follows directly from [15, Theorems 1.1, 2.1 and 2.2]. 


4 Operational matrices 


In this section, the operational matrix of fractional integration in the Riemann— 
Liouville sense and product operational matrix for the SHPs vector ®(t) will 
be derived. To this end, the inner product of the SHP H,,(t) and t* is derived 
explicitly. Then, by using this explicit formulation, the operational matrix of 
fractional integration and product operational matrix can be obtained easily. 


Lemma 1. Suppose that, for a positive real number s, the inner product of 
the nth Hahn polynomial H(t) and t® is denoted by A(n,s). It can be derived 
as 


A(n,s) = 


| 
= 

3 
— 
loa 
Ne 
oH 
w 
Ww 

€ 


henrst*T (atr+1)r(p+N—r+l1) 
Tir+DP(N—-rt+)l(at)rG+l)’ 


r=0 k=0 
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in which hyn is defined in relation (5). 


Proof. From the definition of the discrete inner product (.,.),, in (7), we have 
‘a i a 
A (n, 8) ils =w wo (=) ee 7 2, Qn(r)r w(r 


&; (ror) () GE) 


1 “ a, (at+r+1)C(@+N-r+1) 
WUT (r+ (N-r4+1)P(at+1)P(6+1) 


(11) 


Theorem 2. Let ®(t) be the (N +1) SHPs vector defined in (9). The 
fractional integration of order v for this vector can be expressed as 


I’G(t) ~ P B(t), (12) 


where P™ is an (N +1) x (N +1) matrix and its (i,j)th element can be 
defined as 


as hei N*T (k + 1) A(j,k tv) 


, 47=0,1,...,N. 
T(k+v+1)7r(j) ad 


Proof. The ith element of the vector ®(¢) is H;_1(t) and its fractional integral 
of order v can be derived as 


IY Hy_1(t) = I°Qi-1 (Nt) 


i-1 i-1 
hy,i_-1N*T (k +1) 
Vv k kji-1 v+k 
= (3 init) = as aes a . (13) 


By expanding t**” in terms of SHPs, we get 
N 
pom S- 8,444; (t), (14) 
j=0 
in which 6,,; can be derived by using Lemma 1 as follows: 


1 _ AG, k+v) 
79) ” mG) 
Now, by inserting (14) and (15) in (13), we get 


Bye = — (H(t), t°*”), = (15) 
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, TLS hee NTR ENG, BAY) 
. n= (¥ : Tik+y +l) Hit) 


and this leads to the desired results. 


Theorem 3. Let ®(t) be the (N+1) Hahn polynomials vector defined in (9) 
and let V be an arbitrary (N +1) vector. Then, it is convenient to write 


&(t)®7 (t)V = V(t), (16) 


where V is the (N +1) x (N +1) product operational matrix and its (i, j)th 
element can be defined as 


~ 1 


N 
Visa gti = SOV ORG Oe. 47 =O,1ja1.,0. 
=0 


mH) 4 


Proof. The product of two Hahn polynomial vectors ®(t) and ®7(t) is an 
(N +1) x (N +1) matrix as follows 


Ho(t)Ho(t) Ho(t) Vines ) 
: Hy(t)Ho(t) Hy(t)Hi(t) ... Hy(t)Hy/(t) 
&(t)67 (t) = 


As a result, relation (16) can be rewritten as: 


N N 
S> Ag (t)Hi(t)Views = 5° Ay (t) Vier ets, 1=0,1,...,N. 
k=0 k=0 


Now, by multiplying both sides of the above equation by H,(t) and using the 
defined inner product in (6), we get 


S- (Hy, (t). A; (a), H; (t)) Viet = (H; (t), A; (t)) Viatj+is a, J = 0, 1, sh Re yl 


N 
k=0 


Therefore, the matrix V may be given by 
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5 Description of the proposed method 


In this section, SHPs along with their operational matrix of fractional inte- 
gration will be applied to solve TDFOCPs. Consider the optimal control of 
time-varying linear system (1)—(4) with delays in state and control and with 
quadratic performance. First, we approximate all functions involved in the 
problem as follows: 


D’ x(t) ~ XT O(t), 17) 

u(t) ~ UT @(2), 18) 

r(t) ~ RT O(t), s(t) ~ ST O(t), 19) 

a(t) ~ A &(t), c(t) ~ C7 H(t), 20) 

b(t) ~ BT O(t), e(t) ~ E’ G(t), 21) 

f(t) = FPO), g(t) = GT Ot), m(t) ~ GZ O(), 22) 


where X and U are unknown vectors to be determined, while R, S, A, C, B, 
E, G,, Go, and F are known SHPs coefficient vectors that can be computed as 
explained in (8)-(10). By applying the fractional Riemann—Liouville operator 
I” on both sides of the relation (17), we have 


a(t) ~ IYXT O(t) + g1(0) = X7 P B(t) + GT S(O), (23) 


in which P™ is the operational matrix of fractional integration of the SHPs 
vector ®(t) derived in Theorem 2. Hence, by using (23) and (18), the delay 
state function a(t — 7) and control function u(t — 7) can be derived as 


GT@(t—7), O0<t<r, 
u(t—T)= 
(XT PO +d?) O-—7), r<t<1, 
and 
GZO(t—n), O<t<n, 
u(t—n) = 


UTS(t—n), n<t<l. 


These delay functions can be expanded by the Hahn polynomials as 
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a(t—7) = X7TPMO,6 (t) + VG (8), (24) 

and 
u(t—) =U7QO,,® (t) + V,6 (t), (25) 


where V; and V,, are known (NV + 1) vectors, whereas 0, and 2, are known 
(N +1) x (N +1) matrices, which are depend on the delay parameters 7 
and 7, respectively (please, see details in Appendix A). Now, by substituting 
(17)-(25) in dynamical system (4), we get 


XTO(t) = XT P™ G(t)O7 (t) A + G7O(0) AT O(t) + UT G(t)O7 (t)B 
+X? PMO, B(t)67 (t)C + VI G(t)G7(4)C + UTO, 8(t) 67 (t)E 
+V,P@(t)®7 (t)E + FT O(t). 


By using the product operational matrix as described in Theorem 3, it is 
convenient to rewrite this relation as follows: 


XTO(t) = (xTPOA EGTOO)AT HUT REX POOC 
4+VP64UT0, B+ V2 B+ FT) &(t). 
The above expression is satisfied for all t in the interval [0, 1]; therefore 


X? = XT PMA+ GTO(0)AT +UTB + XPT PMO,E 
4+V/C+U'Q,E+V, E+F". 


Now, from this relation, the unknown vector U can be derived as a function 
of unknown vector X, by solving the following linear system: 


UTA=Y, (26) 
where 
a (B “i OnE) (27) 
and 


Y = X7PMA+ GTO(O)AT + XTPMO,C 
+570 +S) E+ F — X?. (28) 


By inserting the obtained vector U in (18), the unknown control function 
u(t) can be derived as a function of the unknown vector X. Therefore, the 
performance index can be derived as: 


1 
J fxo,21,-..,0n] = f 0 (t, X)dt, (29) 
0 
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where 
6 (t,X) = [(xTP aK) if a(t) ©" ()R+ (UT) OT (L)S}. 


Moreover, by applying the Gauss—Legendre quadrature formula in the in- 


terval [0,1], the performance index J [%o,21,...,2nN] can be approximated 
as 
J (x0, 21,...,0N] ~ >_ we® (te, X), (30) 
k=1 


where wz and t, are Gauss—Legendre quadrature nodes and weights, respec- 
tively. Finally, the necessary conditions for the optimality of the performance 
index are 

OF 


Ox, 
which constitute a system of algebraic equations for unknown vector X. By 
solving this system and determining the vector X, the optimal state function 


x(t) and control function u(t) can be approximated by inserting vector X in 
(18) and (23), respectively. 


0, i=01,...,N, (31) 


Algorithm 1 Algorithm in pseudo-code format 

Inputs: N, a, 6, v,7r, 7, 7, the functions r(t), s(t), f(t), gi (t), g2(t), a(t), c(t), and b(t), e(t). 
Step 1: Construct the Hahn polynomials vector ®(¢) and weight function w(t) from (5) 
and (9). 

Step 2: Define an unknown (N +1) vector X and approximate D” x(t) in (17). 

Step 3: Compute the fractional operational matrix P) form Theorem 2. 

Step 4: Compute the coefficient vectors R, S, A,C, B, E, Gi, G2, and F from (19)—(22). 
Step 5: Compute the product operational matrices A,€ and B, E from Theorem 3. 


Step 6: Compute the matrices 0, and Q,, and vectors V; and V, from Lemma 2 in Ap- 
pendix A. 

Step 7: Compute the matrix A and vector Y from (27) and (28). 

Step 7: Compute the vectors U from (26). 

Step 8: Insert the obtained vectors U in (18) to approximate the control function u(t). 
Step 9: Compute the performance index J in (29). 

Step 10: Apply the Gauss—Legendre quadrature formula to approximate J in (30). 
Step 11: Constitute the system of algebraic equations in (31). 

Step 12: Solve the obtained system of algebraic equations in Step 11 to get the unknown 
vector X. 

Outputs: The approximate state and control functions x(t) ~ X74(t) and u(t) ~ 
UT &(t). 
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5.1 Main features of the proposed method 


A direct numerical method based on discrete orthogonal polynomials is pro- 
posed to approximate solution of time-delay fractional optimal control prob- 
lems. As these polynomials are orthogonal with respect to a discrete norm, 
implementation of the proposed numerical method is more efficient and less 
complex in comparison to similar methods that in which continuous poly- 
nomials are used. In the proposed method, the fractional derivative of an 
unknown state function is approximated by using the discrete polynomials. 
Then, the operational matrix of fractional integration together with the dy- 
namical system is used to approximate the control function as a function of 
state function. Consequently, the need for the direct approximation of con- 
trol function and Lagrange multiplier method are eliminated. The numerical 
algorithm of the proposed method in pseudo-code format is designed in the 
following section. 


6 Illustrative examples 


In this section, some illustrative examples are considered to investigate the 
efficiency and accuracy of the presented SHPs method. In all numerical ex- 
amples, the orthogonal SHPs over the interval [0,1] with a = 6 = 1 are used 
to approximate the solution of (1)—(4). For the problems defined in different 
intervals, first, the interval has been transformed into [0, 1]. All the numerical 
simulation have been done by using Maple 17 with 20 digits precision. 


Example 1. In the first example, we consider the following TDFOCP [20]: 


2 

i= >| [x?(t) + u?(t)] dt 
2 Jo 

subjected to the dynamical system 


D’ a(t) =ay(t)a(t-1)+u(t), O<t<2,0<v<1l, 
a(t)=1, -1<t<0, 


where a,(t) is a known continuous function. For a,(t) = 1, v = 1, and 
N = 12, Table 1 provides a comparison between the obtained optimal value 
J and those derived by Walsh series [33], average approximation [5], Bernstein 
polynomials [37], Boubaker functions [34], shifted Legendre polynomials [7], 
Bernoulli wavelet method [35], and linear programming [20]. From the pre- 
sented results of Table 1, it is clear that the numerical results obtained by 
the presented SHPs method have a good agreement with the true solution 
given by the average approximation [5], Walsh series [33], and linear pro- 
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gramming [20]. Furthermore, it is easy to conclude that the reported values 
in [7,34, 35,37] are not in good agreement with the other numerical results. 
In order to investigate the convergence issue of the proposed algorithm, the 
obtained optimal values of J for different values of N and vy are presented 
in Table 2 and compared with the reported values in [7]. For a1(t) = t, the 
proposed method is employed to solve these TDFOCPs for various values of 
vy. In Table 3, a comparison is made between the obtained optimal values 
of J and the derived values by linear programming method [20]. From the 
presented results in Table 3, we can conclude that the derived optimal val- 
ues J are in a good agreement with the numerical solutions reported in [20]. 
Finally, the graph of approximate state function x(t) and control function 
u(t) with N = 12 are plotted in Figures 1 and 2, for ai(t) = 1 and a,(t) =t, 
respectively. From the obtained results, it is clear that the proposed SHPs 
method is efficient and accurate in solving TDFOCP and that, as the frac- 
tional order v approaches 1, the approximate solution converges to that of 
integer-order problem. Moreover, the presented results in Tables 2 and 3 
indicate that numerical results converge well as the number of basis function 
N increases. 


Table 1: Comparison of the performance index values obtained by different methods for 
ai(t) = 1 and v = 1 (Example 1). 


Method Performance index value J 

Presented method 1.64731019 

Linear programming [20] 1.64886527 
Average approximation method [5] 1.6419 
Walsh series [33] 1.6497 

Legendre polynomials [7] 0.4727464 
Bernstein polynomials [37] 0.6381 

Boubaker functions [34] 0.00002674 
Bernoulli wavelet [35] 0.3048 


Table 2: The performance index J for a1(t) = 1 and different values of v (Example 1). 


vy=0.7 vy=0.8 yv=09 v=0.95 v=0.99 

N=6  1.5581796 1.5964603 1.6266313 1.6384862 1.6464170 
N=8  1.5487500 1.5907820 1.6239023 1.6368166 1.6454020 
N=12 1.5410926 1.5862331 1.6239023 1.6357471 1.6448406 
(7] - 0.4985242 0.5021900 - 0.4778890 
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Table 3: The performance index J for a1(t) = t and different values of vy (Example 1). 


vy =0.7 y=0.8 y=0.9 vy = 0.95 v=1 
N=6_ 0.92467418 0.97005408 1.01321336 1.03283445 1.05078685 
N=8 _ 0.90796836 0.95993974 1.00820712 1.02968701 1.04905258 
N=12 0.90014337 0.95537457 1.00625088 1.02863607 1.04863866 


[20] = 1.08069403 1.06577389 = 1.05137240 
—— o=0.7 —*— a=0.8 — - o=0.9 — o-0.95 —— 0=0.7 —*— o=0.8 — - o=0.9 — o=0.95 
— -cel — 7oFl 


1.04 


u(t) 


0.64 


0.54 


Figure 1: The approximate solutions (¢) and u(t) for a(t) = 1 and N = 12 (Example 


i). 


—— o=0.7 —*— o-0.8 — - o=0.9 —— a-=0.95 —— 0=0.7 —*— o=0.8 — - o=0.9 —— o=0.95 
— -col — 7oFl 


1.0 


u(t) 


Figure 2: The approximate solutions x(t) and u(t) for ai(t) = t and N = 12 (Example 


1). 


Example 2. In this example, we consider the following TDFOCP [34,35] 
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1 


l= 5/ [x?(t) + u?(t)] dt 


subjected to the dynamical system 


D’ x(t) =ta(t)+a(t-1)+u(t), 0<t<2,0<v<1, 
a(t)=0, -1<t<0. 


The SHPs and presented technique in section 5 are used to solve this problem 
for various values of v. For v = 1, Table 4 provides a comparison between 
the results of our proposed method and those derived by FD method [22], 
recursive shooting method [21], line-up competition algorithm [10], Legendre 
multiwavelets [25], Walsh series [33] and Bernoulli wavelet [35]. From Table 
4 we can see that there is good agreement between the obtained results and 
the those reported in Refs. [10, 21,22]. To investigate the convergence issue 
of the proposed algorithm, the obtained optimal values of J for different 
values of N and v are presented in Table 5. Moreover, Figure 3 shows the 
approximate state and control functions x(t) and u(t) for various values of 
vy and N = 12. Based on the obtained results, it is clear that proposed 
SHPs method is efficient and accurate for solving such problems and the 
approximate solutions converge to the exact solution as fractional order v 
approaches integer order 1. Moreover, it is easy to conclude that numerical 
results converge well as the number of basis function N increases. 


Table 4: Comparison of the performance index values obtained by different methods for 
v = 1 (Example 2). 


Method Performance index value J 

Presented method 4.792110320 
FD method [22] 4.79678 
Recursive shooting method [21] 4.79682 
Line-up competition algorithm [10] 4.7976 
Legendre multiwavelets [25] 5.1713 
Iterative dynamic programming [11] 5.0674 
Walsh series [33] 6.0079 
Bernoulli wavelet [35] 2.0481 


Table 5: The performance index J for different values of v (Example 2). 


vy=0.7 vy=0.8 vy =0.9 vy = 0.95 v = 0.99 
N=6 4.04799977 4.32580075 4.58104543 4.69543943 4.77884109 
N=8 4.01176960 4.30195764 4.56856743 4.68746981 4.77382051 
N=12 3.98516242 4.28472972 4.56062884 4.68299652 4.77143293 
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== 0=0.7 — "= 00.8 — ~ 0=0.9 — 0=0.95 == 00.7 —* = o=0.8 — ~ 0=0.9 — 0=0.95 
— _-oFl 
iG 


0.87 


u(t) 


0.24 


Figure 3: The approximate solutions x(t) and u(t) for N = 12 (Example 2). 


Example 3. In this example, we consider the following TDFOCP with 
delays in state and control [7,34, 35] 


D’ x(t) = —a (t)+a(t— 3) +u(t)— su(t- 3), 0<t<1, 0<v<l, 
a(t) =1, “gst, 


The proposed SHPs method has been employed to approximate the solution 
of this TDFOCP for different values of v. For vy = 1, this problem has been 
solved by using the Hybrid functions [17,28] and Bezier curve method [14]. 
Table 6 provides a comparison between the results of our proposed method 
and those reported in [7,14,17,28,34,35]. From Table 6, we can see that there 
is a good agreement between our obtained results and the previously reported 
numerical results in [17,28]. Also, this Table indicates that the reported 
values by using Boubaker polynomials [34], Shifted Legendre polynomials 
[7], Bernoulli wavelet [35], and Bezier curve method [14] are not in a good 
agreement with our proposed method and the valid values given in [17, 28]. 
To investigate the convergence issue of the proposed algorithm, the obtained 
optimal values of J for different values of N and v are presented in Table 7 
and compared with the presented results in [7]. Moreover, Figure 4 shows 
the approximate state and control functions x(t) and u(t) for various values 
of and N = 12. Based on the obtained results, it is clear that the proposed 
SHPs method is efficient and accurate for solving such problems and that 
the approximate solutions converge to the exact solution as fractional order 
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v approaches integer order 1. Moreover, it is easy to conclude that numerical 
results converge well as the number of basis function N increases. 


== 0=0.7 —-— 0=0.8 —~- a=0.9 — a=0.95 == 0=0.7 —-— 0=0.8 —~- a=0.9 — 0=0.95 
— -oFl 


u(t) 


Figure 4: The approximate solutions x(t) and u(t) for N = 12 (Example 3). 


Table 6: Comparison of the performance index values obtained by different methods for 
vy = 1 (Example 3). 


Method Performance index value J 
Presented method method 0.373230152 
Hybrid Legendre functions [28] 0.37311241 
Hybrid Bernoulli functions [17] 0.37310517 
Bernoulli wavelet [35] 0.1027 
Bezier curve method [14] 0.4220497643 
Boubaker polynomials [34] 0.04553 
Legendre polynomials [7] 0.01451 


Table 7: The performance index J for different values of v (Example 3). 


vy=0.7 vy =0.8 y=0.9 y=0.95 v=0.99 
N=8  0.3393198 0.3496035 0.3600181 0.3652627 0.3694686 
N=12 0.3415913 0.3515520 0.3613160 0.3661582 0.3700268 
N=16 0.3549220 0.3549220 0.3637609 0.3683158 0.3720233 
Ref. [7] — 0.0126951 0.0116573 — 0.0139046 


Example 4. Finally, consider the following TDFOCP [7, 14] 
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subjected to the dynamical system 


D’ x(t) =ax(t)t+u(t—)+ul), O<t< 4, O0<v<l, 
u(t) =0, -G<t <0, 
x(0) = 1. 


Basin and Gonzalez [6] and Ghomanjani et al. [14] solved this problem for 
v = 1 by using the Bezier curve and Linear-quadratic method, respectively. 
Here, this TDFOCP has been solved by using the proposed Hahn polynomial 
method for variose values of v. In Table 8, a comparison is made between 
the optimal values of J obtained by the SHPs method and those derived by 
the Bezier curve method [14], shifted Legendre polynomials [7], and Linear- 
quadratic method [6]. From this Table, we can realize that the SHPs method 
is efficient and that the obtained results are in good agreement with the 
presented results in [6,14]. However, the reported results in [7] are not in 
good agreement with other methods. For noninteger values of v, the obtained 
optimal values of J are presented in Table 9 and compared with the value 
of J derived by shifted Legendre polynomials [7]. Finally, Figure 5 shows 
the approximate state function x(t) and the control function u(t) for various 
values of y and N = 12. Based on the obtained results, it is clear that the 
proposed SHPs method is efficient and accurate for solving such problems and 
that the approximate solutions converge to the exact solution as fractional 
order @ approaches integer order 1. Moreover, it is easy to conclude that 
numerical results converge well as the number of basis function N increases. 


Table 8: Comparison of the performance index values obtained by different methods for 
vy = 1 (Example 4). 


Method Performance index value J 
Presented method method 0.1554640279 
Bezier curve method [14] 0.1565866913 
Linear-quadratic [6] 0.1563 
Legendre polynomials [7] 0.0143671 


Table 9: The performance index J for different values of N and v (Example 4). 


vy=0.7 vy =0.8 vy=0.9 y=0.95 v=0.99 
N=6 _ 0.1594198 0.1574747 0.1553910 0.1543186 0.1534526 
N=8 _ 0.1607079 0.1584826 0.1561770 0.1550124 0.1540808 
N=12 0.1616864 0.1592375 0.1567619 0.1555286 0.1545484 
Ref. [7] a 0.0126951 0.0116573 - 0.0139046 
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== 0-0.7 —° = 0=0.8 —~ 0=0.9 — a=0.95 == 0-0.7 =" = 0=0.8 —~ o=0.9 — 00.95 
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Figure 5: The approximate solutions z(t) and u(t) for N = 12 (Example 4). 


7 Conclusion 


A direct numerical method based on the discrete SHPs was proposed to ap- 
proximate solution of TDFOCPs with quadratic performance index. First, 
the operational matrix of fractional integration in the Riemann-—Liouville 
sense and product operational matrix for the shifted Hahn polynomials were 
derived. The fractional derivative of unknown state function was approxi- 
mated by using the SHPs. Then, the operational matrix of fractional inte- 
gration together with the dynamical system were used to approximate the 
control function as a function of state function. Finally, these approximations 
were put in the performance index and necessary conditions for optimality 
transform the under consideration TDFOCPs into an algebraic system. Nu- 
merical results confirm that the presented SHPs method is accurate and 
effective. 
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Appendix A 


Lemma 2. Consider the function 


GT ®(t—7), 0<t<r, 
a(t—T)= 
XTPMO(t—7)+ dT O(t— 7), r<t<1, 


it can be expanded into SHPs as follows: 
a(t—7) = X7PWMO, O(t) + V,B(2), (32) 


in which Q, and S, are, respectively, (N +1) x (N +1) and (N +1) x1 
matrices defined as follows: 


Proof. By expanding the function x(t — 7) expanded by the SHPs, we get 


N 
a(t—7) =) > 4H; (t) = ZT &(t), 
1=0 


where ith element of the coefficient vector Z can be derived as 


N 


wn Mepiie « BSe(E) m(B)e(B). eet 


k=0 


Let J be the greatest integer such that |< 7N. Then 


a= Ey a(S—s)m(B)o(E) +222" 9 o(4 2) n(h) 


Now, let us define 
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3 Ho (1) HH) 0(4) 
1 | Sm (h- D(H) OCF) 
Xi = Ty | k=0 , 
mt) |. 
3 Bw (A 7) i (H) 0(K) 
and 
HoH 7) Hi (H)0(4) 
mia | ge CB) CR) 
a | 
Hw (fe 7) Hi (#) OC) 
k=I4+1 


then, the element z; can be obtained as 
z= Gly, +07, + XP PMT. 
Consequently, we get 
a(t —7) = XTPMO,G(t) + VF H(t), 
where V; is an (N + 1) x 1 vector as 


VE = [GT x0 + d?¥%, GPx t+ d?M%,..., GT xn + dq] 


and Q, is an (N +1) x (N +1) matrix which its ith column is II, that is, 


Ho (4 = 7) Ho (4) 0 (F) - Ho (4-7) Hw (4) 0(4) 
R22 Paz (iy — 7) Ho (w) @ Ge) - ym (4 —7) Hy (%) 0 (4) 
7 ri) | = 


3: Hy (4-7) Ho(#)a(4).. SS Hw (4-7) Hw (4) 0(4) 


k=l+1 k=l41 


Remark 2. In the same way, it can be proved that the 


u(t—) =U7O,® (t) + V,6 (t), 
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where 2, and V,, are (N +1) x (N +1) and (N +1) x 1 matrices, respectively. 


